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The average absorbed dose in glandular tissue is the most appropriate parameter for the assessment of the 
radiation-induced risk during breast imaging. The aims of this work concern: (1) the investigation of the vari- 
ation effect of any related update to photon cross-section data-bases on the computation of the normalized 
glandular dose (DgN) for mammography quality control tests and (2) the proposition of a parameterization 
method leading to provide DgN values function of the breast thickness (T) and the particle energy (E) instead 
of E alone, as normally known. We analyzed the change effect of the photon cross-section data-bases on the 
computation of DgN. Those coefficients, generated using the Geant4 Monte Carlo toolkit, were studied over a 
range of compressed breast thickness of 2-8 cm for monoenergetic (1-120 keV by 1 keV intervals) and polyen- 
ergetic (23-35 kVp by 2kVp intervals) X-ray beams. Moreover, breast tissue composition ranging from about 
0% glandular (about 100% adipose) to 100% glandular (0% adipose) was also covered. The successful param- 
eterization of DgN look-up table function of the breast thickness and energy, will compact its analytical form 
without loss of accuracy. All parameterization fits resulted in r? values of 0.999 or better. 
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I. INTRODUCTION 


Since the glandular breast tissue is regarded as the most 
radio-sensitive tissue, breast dosimetry has been considered 
as an important basis for assessment of radiation risk of pa- 
tients undergoing mammography [1]. The average absorbed 
dose in the glandular breast tissue, known as the mean glan- 
dular dose (MGD), is preferred for radiation risk assessmen- 
t. It can be estimated by measuring the entrance air kerma 
and by using normalized glandular dose (DgN) conversion 
coefficients for a given breast thickness and glandularity [2]. 

Among the most prominent factors influencing DgN val- 
ues was the photon fluence per exposure factor (photon quan- 
ta per mm? per mR). Such factor depends on the mass en- 
ergy absorption coefficient tables for air (en /p) which can 
be extracted from existing data bases, i.e. MCPLIB [3], 
XCOM [4]. Moreover, the use of MCPLIB cross section 
leads to 10% higher conversion coefficient values than the 
use of XCOM data, as pointed out by Zoetelief et al. [2]. 
Thus, the need to update existing monoenergetic DgN(E), 
look-up tables reveals importance for any occurring change 
in the (Hen/p) tables, and the immediately dependent Monte 
Carlo simulation programs used to generate such data. 

On the other hand, due to the large amount of tabulated 
DgN values used for optimizing mammographic procedures, 
provided by many authors [5—10], the requirement of param- 
eterization is obvious in order to simplify direct calculation of 
MGD after measuring the incident air kerma. Thus, we can 
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refer to many works focused on that issue [11-13]. Neverthe- 
less, introduction of breast thickness as a second parameter 
other than energy for parameterization of data-bases does not 
seem well studied, which can lead to more generalizing the 
fitting equations. 

In this paper, we study effect of the photon fluence per 
exposure parameter on the calculation of DgN values; fact 
that seems essential after observation of variation of (Jen/p) 
values between the derived data from Hubbel et al. [14] and 
those used by Boone et al. [11]. The data of monoenergetic 
DgN conversion coefficients are parameterized so that sim- 
ple custom function can be easily handled on a spreadsheet 
to calculate DgN values for a given set of input parameter- 
s (photon energy, breast composition and thickness). Then, 
we adopt a mathematical fitting method, where we include 
the breast thickness as a second parameter and compare their 
accuracy to predict direct computational results. Finally, we 
check ability of the method to expect MGD for clinically rel- 
evant study cases, as carried out by many authors [15], of X- 
ray tube anode/filter combinations of Mo/Mo (30 um), W/Rh 
(50 um) and W/Pd (50 um), 23-35 kVp, breast thickness of 4 
cm and glandularity of 50%. This work can be considered as 
a continuation of many recent works around the MGD calcu- 
lation topic [16-18], performed by medical physicist as part 
of mammography quality control tests. 


Il. MATERIALS AND METHODS 


There are many Monte Carlo simulation packages that 
can be used for estimating the MGD, with different advan- 
tages and disadvantages. Accurate and versatile general- 
purpose simulation packages such as Geant3 [19], EGS4 [20], 
MCNP [21] and most recently Geant4 [22] include well vali- 
dated physics models, geometry-modeling tools and efficient 
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visualization utilities and do not require any additional effort 
to be tailored for medical imaging application. For the current 
study we proceed with the Geant4 toolkit. 

To afford an updated data set of DgN coefficients for mam- 
mography, covering the energy interval of 1-120keV, for 
thicknesses of 2—8 cm and glandularities (0%-100%) of com- 
pressed breast, we considered the same experimental set- 
up used by Boone et al. [11]. The studied energy interval 
gives the possibility of using look-up tables for standard (up 
to 40kVp) and dual-energy (up to 120kVp) mammography 
imaging. We developed a Monte Carlo simulation program 
to generate the data. Mathematical fitting procedure was fol- 
lowed to parameterize the large amount of data. 

According Refs. [1, 23, 24], the MGD calculation is highly 
influenced by spatial distribution of glandular tissues with- 
in the breast. Thus, the current study is dedicated to further 
enhance the quality assurance and quality control procedures 
involving equipment performance, comparisons of X-ray ma- 
chines efficiency and mammography dosimetry protocols. 

For possible larger compressed breast thickness of 8 cm, 
the same parameterization procedure should be followed with 
the penalty of complicating equation and some extra fitting 
parameters will be needed. 


A. Simulation procedure 


In order to derive DgN(E) values, we used the photon quan- 
ta per mm? per mR parameter, K(E), for a given energy E (in 
keV) of incident photon beam in a medium of air as derived 
by Johns et al. [25] 


K(E) = 54300/[E (Hen/)air] o) 


where, (Hen/P)air, in cm?/g, is the mass energy absorption 
coefficient in air. 
The DgN(E) coefficient, in mGy/mGy, is given by [26, 27] 


-17 
DeN(E) = 1.8352 x 10 Ea G{E) KE) (2) 
fep(T — 2Tskin)(R2/ Ri)? 

where, Edep is the energy deposited in the breast tissue, fg 
is the glandular fraction of the breast, p is the breast density, 
R, = 8cm is breast radius, Rə = 7.6cm is breast tissue ra- 
dius excluding the skin, Tyxin is the skin thickness, T is the 
breast thickness and G(E) is a factor describing the propor- 
tion of the absorbed dose of glandular tissue to the overall 
breast tissues and given by 


= fg (Hen/P)g 
lel Gs), =F ua 


where, (Hen/P)g and ({ten/)a are the mass energy absorption 
coefficients for glandular and adipose tissue, respectively, for 
the given energy. The unit of the constant 1.8352 x 10717 is 
(mGy g)/keV. 

Using the recent mass energy attenuation coefficients of 
photons through air, provided by NIST laboratory [14], a 
Monte Carlo-based simulation program was tailored for cal- 
culating absorbed energy in the breast model for a given 
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configuration of breast tissue composition, thickness and for 
monoenergetic photon beam. Then, an additional procedure 
was added to that program, in order to compute the G(E) fac- 
tor during the simulation run-time, which in turn outputs the 
DgN(E) coefficients for each studied case. 

The Geant4-based simulation program was tailored to 
model monochromatic X-ray emission from a “point-like” 
source and to follow particle transportation through the 
breast. The physical processes used were the photoelectric 
effect, the Compton and the Rayleigh scattering for photons 
and the bremsstrahlung, the ionization and the “multiple scat- 
tering” for electrons. The production thresholds were set to 
1 keV for photons and 10 keV for electrons. From the existing 
suite of physical packages, ready to use within the simulation 
program, the “PhysListEmStandard” was selected. Elemen- 
tal compositions of adipose, glandular and skin tissue were 
based on the work of Ref. [27]. 

For each set of (monoenergetic) X-ray energy, breast size 
and breast composition, a simulation run was executed using 
10° entrant photons providing sufficient statistical precision. 
In order to simplify analytical form of the fitting functions 
to be used and due to the small contribution to MGD param- 
eters, we omitted the interval of energy of 1-11 keV during 
the fitting stage of the database. Notice that such contribution 
does not exceed 1% of the MGD for the studied X-ray spectra. 
DgN(E) coefficients for 11-120 keV X-rays were computed. 
At least a total number of 1582 (113 x 7 x 2) runs were carried 
out for the data-base generation purpose covering 113 ener- 
gy values (11-120 keV), and seven thicknesses (2-8 cm) for 
glandularities of 100% and 50%. The absorbed dose by the 
breast tissue was recorded for direct derivation of the DgN(E) 
coefficients (in the SI units of mGy/mGy). 

As we followed the same simulation setup used by Boone 
et al. [11], the same source-to-detector distance of 60 cm was 
used. Changing the distance to 80cm or 40cm can deviate 
the beam energy at the vicinity of the breast by less than 1%, 
which is logical due to the fact that the addition or subtraction 
of a 20cm layer of air medium can affect the primary beam 
trajectory by 0.4% and 0.06% for 30 keV and 100 keV pho- 
ton energy, respectively. Moreover, the assumption of model- 
ing the geometry of the X-ray tubes focal spot as a point-like 
rather than a disc is valid to about 1%. The replacement of a 
point-like source with a disc of 40 um or 300 um (as found 
in clinical realistic tests), deviate the DgN(E) with 0.213% or 
0.444% for a given couple of photon energy, glandularity and 
thickness of (50 keV, 50%, 2cm) or (10 keV, 100%, 8 cm), 
respectively. Whereas the focal spot size affects contrast of 
the images [28]. 

We can consider the homogeneous distribution of the glan- 
dular tissues within the breast as a limitation of the current 
model, as mentioned by Dance et al. [29] that the absorbed 
energy by the glandular tissue depends on the glandular tis- 
sue position in the breast. 

The simulation was carried out using Geant4 (V9.4.p01) 
under a Red Hat Enterprise Linux 5 workstation on an Intel 
Core i7 (CPU) running with a 4 GB RAM at 3.40 GHz. The s- 
tatistical uncertainty (2c) associated with all the Monte Carlo 
calculations presented in this work is less than 1%. 
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B. Parameterization procedure 


The following parameterization procedure was followed to 
reduce the large amount of DgN(E) values while retaining the 
prediction precision to an acceptable level of accuracy and us- 


DgN(E) 


a+cln(T) + eln(£) + gln?(T) +iln?(£) + kn(T) In(£) 
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ing simple analytical form of the fitting equation. For accura- 
cy of the proposed method, the r? statistical coefficient value 
was calculated. We fitted the overall data set of DgN(T,E) 
using Eq. (4) for the energy interval of 11-120 keV, for glan- 
dularities of 100% and 50% separately 


(4) 


And the MGD was calculated using Eq. (5): 


opm K(E) x DgN(E) x (E) x ¢(E) 


MGD = É, 
Deas SA OE) 


» ©) 


where, (E) is the photon fluence per exposure in units of 
mR per photons per mm? and Emin and Emax are the minimum 
and the maximum energy spectrum limits. 

Physically spoken, it is clear from Eq. (4) that any DgN 
coefficient explicitly depends on the coming photon energy 
and the breast phantom thickness [30]. Thus, we prefer to 
fit those values with a bi-variant polynomial depending on E 
and T [13]. The quotient form can be explained by the fact 
that the DgN coefficient is a fraction of the absorbed dose 
by the breast model to the air medium, and each of them can 
be a polynomial energy dependent function. For a relative 
difference between computed and fitted values that does not 
exceed 3% on the MGD calculation for the majority of X- 
ray beam spectra. We were limited to adopt the second order 
polynomial fitting Eq. (4). More precision needs extension 
to higher order, with the penalty of complicating the fitting 
expression form. 

To our knowledge, fitting methods using the parameter E 
was treated similarly by many authors [7, 26, 30], whereas the 
introduction of the parameter T as a second variable within its 
actual form is proposed for the first time herein. 


Ill. RESULTS AND DISCUSSION 


We started with studying effects of photon fluence per ex- 
posure coefficient which, in turn, caused the modification in 
DgN value between those of Boone et al. [11] and those p- 
resented during this work. Moreover, the fitting parameters 
were tabulated for any easy handling of the computational 
model. 

Based on the parameterization procedure, we studied the 
efficiency to reproduce directly simulated MGD for some 
specific cases. We terminated this work with analyzing the 
interpolation and the extrapolation capabilities of such pro- 
cedure. Further work can be conducted towards the use of 
heterogeneous breast tissue composition and within a prone 
breast situation. 


~ 14 bIn(T) + dln(E) + fln?(T) +h? (E) + jIn(T) uE) 


A. X-ray Quanta per mm” per mR computation 


We calculated the photon fluence per exposure conversion 
coefficient, using Eq. (1), for photon energy of 1-120 keV. 
Figure 1 shows the relative difference between K(E) calculat- 
ed based on the most recent data of mass energy absorption 
coefficient of air, provided by NIST [14], and those tabulat- 
ed by Boone et al. [26]. The major difference seen concerns 
photons with energy lower than 75 keV and can reach 12% for 
27 keV, as seen in the inset. For photon energies over 75 keV, 
such difference does not exceed 1%. 


2 Boone 97 
+ NIST $ 3 


“9 20 40 60 80 100 120 
Photon Energy (keV) 


Fig. 1. Comparison of the photon quanta per mm? per mR parameter 
calculated in this work (denoted by NIST) against Boone et al. [26] 
data (denoted by Boone 97). The inset refers to their relative devia- 
tion as function of the energy. 


B. K(E) variation effect on Boone’s look-up table 


Table 1 illustrates a comparison of MGD calculations be- 
tween literature (Boone ef al. [5]) and present study using 
the updated K(E) values. The comparison concerns three 
different anode/filter combinations of Mo/Mo (30 um), W/Rh 
(50 um) and W/Pd (50 um) for 23-35 kVp and HVL of 0.269- 
0.597mm of Al. The breast thickness was 4cm and the 
glandularity was 50%. We observed a maximum relative 
difference between literature and Geant4 data of about 1.9, 
5.0 and 4.2%, and between calculated (parameterized) and 
Geant4 data of about 2.6, 3.8 and 4.1%, respectively, for the 
above anode/filter combinations. 
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TABLE 1. MGD (mGy/mGy) calculated for different combinations of anode/filter, kVp and HVL: using Boone’s method (denoted by Boone 
et al. [5]), directly simulated using Geant4 (denoted by G4) and using the parameterization method (denoted by Calc.). Deviation between 


Boone et al./Calc. and G4 data were tabulated in % 


Anode/filter E AVL 
(kVp) (mm Al) Boone et al. [5] 
Mo/Mo30 23 0.269 0.158 
25 0.295 0.179 
27 0.318 0.195 
29 0.338 0.210 
31 0.356 0.222 
33 0.372 0.233 
35 0.386 0.243 
W/Rh50 23 0.420 0.258 
25 0.462 0.289 
27 0.489 0.307 
29 0.509 0.321 
31 0.527 0.332 
33 0.544 0.344 
35 0.560 0.355 
W/Pd50 23 0.424 0.260 
25 0.478 0.300 
27 0.514 0.324 
29 0.539 0.340 
31 0.560 0.353 
33 0.579 0.365 
35 0.597 0.377 
30 Tw peak —-- KE) 
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Fig. 2. Relative difference curves of DgN(E), K(£) and their sum 
vs. E, between those computed in this work to those derived using 
old cross section data-base and Monte Carlo simulation code found 
by Boone et al. [5, 22]. The two energy intervals (grey bands) in- 
dicate the localization of the main peak of the tungsten (W) and the 
molybdenum (Mo) experimental spectra. 


Referring to the Mo and W experimental spectra in [31], 
we observed that the two characteristic peaks were located 
around 17 keV and 11 keV, respectively, as shown in Fig. 2. 
The spectra covered 23-35kVp. Also, Fig. 2 illustrates the 
curves of the relative error (%) of K(E), DgN(E) and their 
sum. According to Eq. (5) and Fig. 2, we are intended to find 
the MGD relative difference corresponding to the tungsten an- 
ode twice or higher than those for molybdenum. This can be 
explained by coincidence of the tungsten peak and the max- 
imum relative deviation. Furthermore, we can guess that the 
relative deviation corresponding to the rhodium anode spectra 


MGD Deviation (%) 
G4 Calc. Boone et al. [5] Calc. 
0.155 0.151 1.935 2.581 
0.178 0.175 0.562 1.685 
0.197 0.195 1.015 1.015 
0.212 0.212 0.943 0.000 
0.225 0.225 1.333 0.000 
0.236 0.237 1.271 0.424 
0.244 0.246 0.410 0.820 
0.265 0.266 2.642 0.377 
0.299 0.302 3.344 1.003 
0.321 0.326 4.361 1.558 
0.338 0.349 5.030 3.254 
0.349 0.358 4.871 2.579 
0.358 0.371 3.911 3.631 
0.370 0.384 4.054 3.784 
0.267 0.269 2.622 0.749 
0.309 0.312 2.913 0.971 
0.337 0.341 3.858 1.187 
0.355 0.367 4.225 3.380 
0.366 0.377 3.552 3.005 
0.377 0.392 3.183 3.979 
0.389 0.405 3.085 4.113 


TABLE 2. Computer-fit parameters used to derive DgN(E,T), in m- 
Gy/mGy, for Æ = 11-120 keV and T = 2-8 cm. f: glandularity, 
p: (a, .., k). The statistical p-value (r?) was given 


p fe = 100% fe = 50% 

a —0.1463 —0.1811 
b 0.01196 0.006104 
c 0.1547 0.1986 
d —0.5080 —0.5147 
e —0.03587 —0.04615 
f 0.03782 0.04211 
g 0.03284 0.0359 
h 0.09141 0.09939 
i 0.05131 0.06342 
j —0.05943 —0.06876 
k —0.1040 —0.1269 
r? 0.9999 0.9999 


will be lower than those for the molybdenum, for the studied 
kVp interval. 


C. Data set generation and parameterization 


The overall normalized glandular dose data were generated 
using Eqs. (2) and (3) and the developed simulation code, for 
breast glandularities of 100% and 50%. The data set corre- 
sponding to f, = 50% is shown in Fig. 3. Plots correspond- 
ing to fẹ = 100% are similar to Fig. 3. Table 2 illustrates the 
parameters resulting from fitting the overall DgN values and 
using Eq. (4) for energy interval of 11-120 keV. There is no 
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TABLE 3. DgN(E) coefficients, in mGy/mGy, directly simulated (denoted by G4) and parameterized (denoted by Calc) for different glandu- 
larity values of 20%, 40%, 60%, 80% and 90 %. E: energy in keV, T: thickness in cm 


E(keV) T(cm) fa = 20% fa = 40% 
G4 Calc. G4 Calc. 
15 2 0.34 0.34 0.32 0.32 
50 2 1.20 1.19 1.21 1.20 
15 4 0.18 0.18 0.17 0.17 
50 4 1.18 1.19 1.19 1.19 
15 8 0.09 0.09 0.08 0.08 
50 8 1.03 1.03 1.03 1.02 


DgN(mGy/mGy) 


Fig. 3. (Color online) Simulated DgN (mGy/mGy) versus energy 
(keV) and thickness (cm) for 50% glandularity, generated using the 
Geant4-based simulation program. 


need to subdivide the energy and thickness intervals, for pa- 
rameterization purposes. A statistical test p-value of 0.9999 
was seen for both cases. 


D. MGD prediction ability 


In order to investigate the possibility of MGD prediction 
using the current data set of DgN following the fitting method, 
we used Eq. (5) and the photon fluence per exposure spectrum 
corresponding to the anode/filter combinations of: Mo/Mo 
(30 um), W/Rh (50 um) and W/Pd (50 um) generated accord- 
ing to the MASMIP and the TASMIP Boone’s codes [26] 
for different kVp and HVL X-ray tube parameters. The heel 
effect associated to the spectra was not taken into account for 
this study, similar to [16]. Thus the beam is assumed to have 
a uniform intensity and uniform spectral quality within the 
field covering the entire breast. If we took it into account, the 
MGD would vary as it is proportional to the fluence spectra. 

From the data at 35kVp in Table 1, we can see that MGD 
values become close to each other. This is due to the small 
contribution of DgN for low photon energy and high kVp. 
The results are in acceptable agreement with those directly 
simulated. 


fe = 60% fe = 80% fe = 90% 
G4 Calc. G4 Calc. G4 Calc. 
0.30 0.30 0.28 0.29 0.28 0.28 
1.21 1.20 1.22 1.21 1.23 1.21 
0.15 0.15 0.14 0.14 0.13 0.13 
1.18 1.19 1.19 1.19 1.20 1.19 
0.07 0.07 0.07 0.07 0.06 0.06 
1.02 1.02 1.02 1.01 1.02 1.01 


E. Interpolation and extrapolation capability 


The DgN data values in Table 3 demonstrate effectiveness 
of the interpolation and extrapolation procedures, in terms of 
glandularity. The direct calculation using the linear interpola- 
tion and extrapolation technique and the Monte Carlo simula- 
tion were conducted for different cases of glandularity and the 
data are tabulated. We find an acceptable agreement between 
the two sets of data allowing the straightforward interpolation 
and extrapolation of the MGD for any given glandularity from 
about 0% to 100%. 


IV. CONCLUSION 


We calculated the photon fluence per exposure conversion 
coefficient, using Eq. (1), for the interval of photon ener- 
gy ranging from 1 keV to 120keV. A computer simulation 
program was developed to estimate the radiation dose to the 
breast phantom during quality control and quality assurance 
Mammography tests. Monoenergetic DgN coefficients were 
computed for various energies and thicknesses and the data 
have been parameterized and presented in the form of a gen- 
eral equation and parameters depending on the two variables 
of E and T. The Geant4-based computed DgN coefficients p- 
resented in this work provide an update version of those main- 
tained during mammographic dosimetry quality control tests, 
nowadays. The ultimate need to recompute such coefficients 
is related to any update occurring to the mass-energy ab- 
sorption data-base of materials, which in turn influences any 
Monte Carlo simulation outputs. Even though that the relative 
difference found during the calculation of the K(E) parame- 
ter reached 12% for 27 keV, the relative difference occurring 
to the MGD values attained the 5% for the studied cases of 
Molybdenum and Tungsten anode spectra. The parameteri- 
zation procedure allowed us to correctly compute the MGD 
values for any X-ray spectra. Using the linear interpolation 
and extrapolation techniques, we are able to vary the kVp up 
to 120 keV, the breast phantom thickness from 2cm to 8cm 
and the glandularity proportion from about 0% to 100% and 
correctly reproduce MGD values derived from direct simu- 
lation. Such parameterization procedure can be straightfor- 
wardly applied to modern techniques such as breast CT and 
digital breast tomosynthesis, despite compressed breast CT in 
concern. 
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